<?php
function peptide_blast_pre_process($seq, $eval, $org, $ref) {
	$wd = './blast/';
	$seq = blastseqpack($seq);
	openDatabase();
	$blastout = run_blast_nr_peptide($wd,$seq,$eval,$org,$ref);
	$hits = blastparse($blastout);
	closeDatabase();
	
	unlink($seqfile); //delete seqeunce file
	unlink($blastout);
	//$ip = getRealIpAddr();
	//$blastres = str_replace('/', '', $blastout).'.'.$ip.'.'.time();
	//serializeToDb($blastres,$hits);
	
	return $hits;
}

function protein_blast_pre_process($seq, $eval, $org, $ref) {
	$wd = './blast/';
	$seq = blastseqpack($seq);
	openDatabase();
	$blastout = run_blast_protein($wd,$seq,$eval,$org,$ref);
	//parse the blast result to get gi number and evalue
	$hits = array();
	$lines = file($blastout);
	while (current($lines)) {
		$line = current($lines);
		if (strstr($line,"***** No hits found ******")) {
			break;
		}
		if (strstr($line,"Sequences producing significant alignments:")) {
			next($lines);
			next($lines);
			while (current($lines)) {
					$line = trim(current($lines));
					if (!empty($line)) { // hits
							$id = substr($line,0,strpos($line," "));
							$evalue = substr($line,strrpos($line," "));
							if ($hits[$id]==null) { // if already exists
									$hits[$id] = $evalue;
							}
					} else {
							break;
					}
					next($lines);
			}
			break;
		}
		next($lines);
	} // end parse
	closeDatabase();
	
	unlink($seqfile); //delete seqeunce file
	unlink($blastout);
	//$ip = getRealIpAddr();
	//$blastres = str_replace('/', '', $blastout).'.'.$ip.'.'.time();
	//serializeToDb($blastres,$hits);
	
	return $hits;
}

function run_blast_nr_peptide($wd,$seq,$eval,$org,$ref) {
    $v = getP3DBVersion();
    $fasta = $wd."data/nrpep.$org.$ref.$v";
    if (!file_exists($fasta)) {
        write_fasta_nr_peptide($org,$ref,$fasta);
        run_formatdb($wd, $fasta);
    }
    return run_blast($wd,$seq,$eval,$fasta);
}

function run_blast_protein($wd,$seq,$eval,$org,$ref) {
    $v = getP3DBVersion();
    $fasta = $wd."data/pro.$org.$ref.$v";
    if (!file_exists($fasta)) {
        write_fasta_protein($org,$ref,$fasta);
        run_formatdb($wd, $fasta);
    }
    return run_blast($wd,$seq,$eval,$fasta);
}

function write_fasta_nr_peptide($org,$ref,$dir) {
    $query = "SELECT DISTINCT phosphoSequence,sequence,protein,location FROM peptide WHERE 1";
    if (!empty($org)) {
        $query .= " AND '$org'=(SELECT organism FROM protein WHERE protein.id=peptide.protein)";
    }
    if (!empty($ref)) {
        $query .= " AND '$ref'=(SELECT dataSource FROM mass WHERE mass.id=peptide.mass)";
    }

    $fhandle = fopen($dir, "w");
	$result = runQuery($query);
    foreach ($result as $row) {
        fwrite($fhandle, ">{$row['phosphoSequence']}|{$row['protein']}|{$row['location']}\n{$row['sequence']}\n");
    }
    fclose($fhandle);
}

function write_fasta_protein($org,$ref,$dir) {
    $query  = "SELECT id,sequence FROM protein WHERE 1";
    if (!empty($org)) {
        $query .= " AND organism='$org'";
    }

    if (!empty($ref)) {
        $query .= " AND EXISTS (".
                  "SELECT * FROM proteinDataSourceRelation ".
                  "WHERE protein.id=proteinDataSourceRelation.protein AND dataSource=$ref)";
    }
    
    $fhandle = fopen($dir, "w");
	$result = runQuery($query);
    foreach ($result as $row) {
        fwrite($fhandle, ">{$row['id']}\n{$row['sequence']}\n");
    }
    fclose($fhandle);
}

function run_blast($wd,$seq,$eval,$fasta) {
	$seqfile = tempnam($wd."tmp", "seq_");
	$fhandle = fopen($seqfile, "w");
	fwrite($fhandle, $seq);
	fclose($fhandle);

	$blastout = tempnam($wd."tmp", "blast_result_");
	
	$command = $wd."blastall -p blastp -F F -v 100 -b 100 -i $seqfile -o $blastout -e $eval -d $fasta";
	exec($command);
	
	return $blastout;
}

function run_formatdb($wd,$dir) {
    $command = $wd."formatdb -i $dir";
    exec($command);
}

// parse blast result
function blastparse($blastout) {
	$hitCount = 0;
    $lines = file($blastout);
    while (current($lines)) {
        next($lines);
        $line = trim(current($lines));
        if (strstr(line,"***** No hits found ******")) {
            return null;
        }
        if (strstr($line,"Sequences producing significant alignments:")) {
            next($lines);
            next($lines);
            while (current($lines)) { // skip hit list
                if (trim(current($lines))=="") {
                    break;
                }
                next($lines);
            }
            break;
        }
    }
    // hit data begins
    $hits = array();
    while (current($lines)) {
        while (true) {
            $line = current($lines);
            if (!$line) { // end of file
                return $hits;
            }
            if (trim($line)!="") {
                if (strstr($line,">")) {
                    break;
                }
            }
            next($lines);
        }
        // one hit
        $hit = array();
        $line = trim(current($lines)); //>dds#pps|1|10
        $strs = explode('|', substr($line,1));
        $hit["nrseq"] = $strs[0];
        $hit["protein"] = $strs[1];
        $hit["location"] = $strs[2];
        
        next($lines);
        next($lines);
        next($lines);
        
        $line = trim(current($lines)); // Score = 54.3 bits (129), Expect = 3e-010
        $s = strpos($line,"Score = ");
        $line = trim(substr($line,$s+8));
        $s = strpos($line," bits (");
        $hit["score"] = substr($line, 0, $s);
        $line = trim(substr($line,$s+7));
        $s = strpos($line,"), Expect = ");
        $hit["bits"] = trim(substr($line, 0, $s));
        $hit["evalue"] = floatval(substr($line,$s+12));
        
        next($lines);
        $line = trim(current($lines)); // Identities = 10/22 (45%), Positives = 13/22 (59%), Gaps = 3/22 (13%)
        $arr = explode(")",$line);
        $hit["ide"] = substr($arr[0],strpos($arr[0],"(")+1);
        $hit["pos"] = substr($arr[1],strpos($arr[1],"(")+1);
        if (count($arr==2)) {
            $hit["gap"] = "0%";
        } else {
            $hit["gap"] = substr($arr[2],strpos($arr[2],"(")+1);
        }
        
        next($lines);
        next($lines);
        
        $line = trim(current($lines)); //Query: 184 HSLHSYFLLVGD 195
        $s = strpos($line,"Query:");
        $line = trim(substr($line,$s+6));
        $s = strpos($line," ");
        $hit["querystart"] = intval(substr($line,0,$s))-1;
        $line = trim(substr($line, $s));
        $s = strrpos($line," ");
        $hit["queryseq"] = trim(substr($line,0,$s));
        $hit["queryend"] = intval(trim(substr($line,$s)))-1;
        
        next($lines); 
        $line = trim(current($lines)); //T  D+P+FG+  G   GQ  GQ
        $hit["matchseq"] = $line;
		
        next($lines);
        $line = trim(current($lines));
        $s = strpos($line,"Sbjct:");
        $line = trim(substr($line,$s+6));
        $s = strpos($line," ");
        $hit["sbjctstart"] = intval(substr($line,0,$s))-1;
        $line = trim(substr($line, $s));
        $s = strrpos($line," ");
        $hit["sbjctseq"] = trim(substr($line,0,$s));
        $hit["sbjctend"] = intval(trim(substr($line,$s)))-1;
		
        // index of site are ungapped
        $hit["sites"] = getsbjctsites($hit["nrseq"], $hit["protein"], $hit["location"], $hit["sbjctstart"], $hit["sbjctend"]);
		
        if (count($hit["sites"])>0) {
            $hits[$hit['protein'].'.'.$hitCount] = $hit;
			$hitCount++;
        }
        next($lines);
    }
}

function getsbjctsites($nrseq, $protein, $location, $sbjctstart, $sbjctend) {
	$locInPep = getSiteFromNrPepSeq($nrseq);
	$sites = array();
	foreach ($locInPep as $loc) {
		if($loc<$sbjctstart||$loc>=$sbjctend)
			continue;
		$siteLoc = $loc + $location;
		$query = "SELECT id FROM site WHERE protein=$protein AND locationInProtein=$siteLoc AND ptm='p'";
		$result = runQuery($query);
		foreach($result as $row) {
			$sites[$loc-$sbjctstart] = $row['id'];
		}
	}
        return $sites;
}

function getSiteFromNrPepSeq($nrpseq) {
    $curr = strpos($nrpseq, '#');
    while ($curr) {
        $sites[] = $curr - count($sites) - 1;
        $curr = strpos($nrpseq, '#', $curr+1);
    }
    return $sites;
}

// remove non-alphabetical characters
function blastseqpack($seq) {
    $seq = trim($seq);
    if (strlen($seq)==0) return $seq;
    
    $lines = explode("\n",$seq);
    $line = trim(current($lines));
    while (empty($line)) {
		next($lines);
		$line = trim(current($lines));
    }
    if (substr($line,0,1)==">") {
        next($lines);
        $line =trim(current($lines));
    }
    
    do {
        $return_this = $return_this.eregi_replace("[^a-z]","",$line);
        next($lines);
        $line = trim(current($lines));
    } while ($line && !empty($line) && substr($line,0,1)!=">");

    $return_this = strtoupper($return_this);
    return $return_this;
}
?>